Introduction 


PVLD 


o 

fa  M 

§:S»  e  « 


T77" 


G 

o 

*0  *ri 

<D  «-> 


I  W  <  §  -H 
H  O  V, 

w  o  c  2 
M  M  C3  M 
£  H  C  3 
*  «  P  fa 


,  'O  „ 

i  ^  ^ 

|  T3  r-» 

V  S-g  ® 
O  .0 

W  r-f  Cl 

♦»  ’  3  *2  ,$• 

3  .c  g  w 
J3  ,  ffl  ► 

*H  ;  r-l  _ 

**l 


ffl 

,.  Ml  > 
S  *H  <c 

«  <v 


♦> 

w 


ad 


;The  present  report  investigates  different  numerical  methods  of 
solution  to  compute  wave  forces  on  submerged  shell  type  structures  of 
arbitrary  shape.  Although  we  only  dealt  here  with  two  dimensional 
applications,  future  applications  will  involve  three  dimensional  analysis. 
Hence  it  is  of  primary  importance  that  the  technique  under  consideration 
can  be  extended  to  three  dimensional  applications. 


The  report  considers  for  comparison  purpose  the  application  of  the 
finite  element  method  [l]  but  it  is  evident  that  the  technique  can  not 
be  easily  extended  to  general  three  dimensional  problems.  Thus  we  con¬ 
centrate  in  the  boundary  element  method  [2 | [3]  which  is  newly  developed 
technique. 


The  finite  element  method  has  attracted  the  attention  of  the  analysts 
to  largely  due  to  its  property  of  dividing  the  continuum  into  a  series  of 
elements,  which  can  be  associated  with  physical  parts.  The  existing 
literature  or  finite  elements  is  by  now  very  extensive  and  encompasses 
structural  (4]  as  well  as  fluid  flow  jj  j  and  other  types  of  problems.  The 
method  can  sometimes  be  based  on  variational  principles  or  more  generally 
on  weighted  residual  expressions.  Integral  equations  techniques  on  the 
other  hand  were  until  recently  considered  to  be  a  different  type  of 
analytical  method,  somewhat  unrelated  to  other  approximate  techniques  such 
as  finite  elements.  They  become  popular  in  Europe  through  the  work  of  a 
series  of  Russian  authors  such  as  Muskhelishvili  [5J ,  Mikhlin  [6 J , 

Kupradze  [7j;  Smirnov  [8]  but  were  not  very  popular  with  engineers.  A 
predecessor  of  some  of  this  work  was  Kellog  [V]  who  applied  integral 
equations  for  the  solution  of  Laplace’s  type  problems.  Integral  equations 
techniques  were  mainly  used  in  fluid  mechanics  and  general  potential 
problems  and  known  as  the  ’source’  method,  which  is  now  called  an  ’indirect* 
method,  i.e.  the  unknowns  are  not  the  physical  variables  of  the  problem. 

Work  on  this  method  continued  throughout  the  sixties  and  seventies  in 
the  pionerring  work  of  Jaswon  [joj,  and  Symm  (j  Q ,  Massonet  [_\2j ,  Hess  [j  3] 
and  many  others. 

It  is  difficult  to  point  out  precisely  who  was  the  first  one  to  propose 
the  ’direct’  method  of  analysis.  It  is  found  in  a  different  form,  in 
Kupradze *s  book  [?}.  It  seems  fair  however  from  the  engineering  point  of 
view  to  consider  that  the  method  originated  with  the  work  of  Cruse  and 


Rizzo  [_I4j  in  elastostatics .  The  inportance  of  the  direct  technique  is 
not  only  that  it  allows  us  to  solve  the  problem  in  terms  of  physical 
variables,  but  also  that  it  was  the  first  step  towards  a  better  under¬ 
standing  of  the  technique  and  its  relationship  to  other  approximate 
methods,  in  particular  mixed  finite  element  models,  such  as  those  proposed 
by  Pain  [_15^ . 

Since  the  early  1960's  a  small  research  group  at  Southampton  University 
started  working  <>n  the  application  of  integral  equations  to  solve  engineering 
problems.  Hadid's  thesis  published  in  1964  [l6"[  dealt  with  the  applications 
of  integral  equations  in  shell  analysis.  Unfortunately  the  presentation  of 
the  problem,  the  difficulty  of  defining  the  appropriate  Green's  functions 
and  the  parallel  emergence  of  the  finite  element  method  all  contributed  to 
minimize  the  importance  of  this  work.  It  was  not  until  1973  that  the  group 
produced  another  significant  thesis,  i.e.  the  one  written  by  Watson  (j 7j 
and  dealing  with  three  dimensional  elastostatics  problems.  By  then  recent 
developments  in  finite  eleipents  had  started  to  find  their  way  into  the 
formulation  of  boundary  integral  equations.  The  idea  of  using  general 
curved  elements  originated  in  another  Southampton  University  thesis 
presented  by  Lachat  in  1975.  This  thesis  whose  topic  was  suggested 

and  supervised  by  C.  Brebbia  marked  emergence  of  the  boundary  element 
method.  Still  the  question  of  how  to  effectively  relate  the  boundary 
integral  equations  to  other  approximate  techniques  was  still  unresolved. 

This  was  done  by  Brebbia  who  published  a  series  of  papers  on  the  relationship 
of  different  approximate  methods,  eliminating  in  1978  with  the  first  book 
[dj  for  which  the  title  Boundary  Elements  was  used.  More  recently  this  work 
has  been  expanded  to  encompass  time  dependent  and  non-linear  problems  and 
this  gave  origin  to  a  more  advanced  second  boundary  elements  book  [3] . 

Two  inportant  International  Conferences  were  held  at  Southampton  University 
in  1978  and  1980.  The  edited  proceedings  of  these  conferences  -  so  far 
the  only  ones  on  the  topic  -  are  now  standard  references  [l8]£l9^]. 

Although  boundary  elements  is  a  very  powerful  technique  it  is  convenient 
in  applications  such  as  fluid-shell  interaction  problems  to  combine  it 
with  finite  elements.  This  combination  should  be  achieved  satisfying 
fully  compatibility  and  equilibrium  at  the  interfaces  between  fluid  and 
solid,  which  requires  using  the  same  types  of  elements  for  both  solutions. 

The  coupling  may  be  achieved  in  either  of  two  ways;  i)  by  considering 
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the  whole  problem  using  an  equivalent  BEM  approach;  ii)  by  converting 
the  BEM  regions  into  equivalent  FEM.  The  two  approaches  are  described 
in  detail  in  reference  [20J  .  As  the  FEM  is  very  well  established,  the 
consideration  of  the  BE  subregion  as  an  additional  equivalent  finite 
element  matrix  seems  most  attractive.  The  formulation  of  this  'equivalent' 
matrix  used  to  model  the  BE  presents  however,  certain  problems;  at 
sharp  geometric  discontinuities,  there  are  also  discontinuities  of  surface 
tractions  which  require  special  attention  and  the  equivalent  FE  matrix 
formed  is  not  inherently  syninetric,  unlike  the  classical  FE  approach. 

A  technique  which  overcomes  these  problems  and  provides  an  acceptable  FE 
type  formulation  using  the  BEM  method  has  recently  been  presented  by 
Georgiou  and  Brebbia.  An  early  paper  involving  a  crude  way  of  symmetrizing 
the  matrices  has  been  published  by  Zienkiewicz  [22]]  et  al,  but  this 
work  has  been  largely  superceded  by  a  recent  paper  by  Mustoe  et  al  [23], 
who  presented  an  interesting  way  of  overcoming  some  of  the  difficulties 
due  to  the  non-symmetric  way  in  which  the  integrations  are  carried  out. 

It  involves  weighting  the  standard  integral  equation  relationships  with 
the  actual  shape  functions  used  for  the  boundary  variables.  This  in 
effect  means  that  instead  of  applying  a  point  source  to  form  each  equation 
a  distributed  source  on  either  side  of  the  node  is  applied  and  a  numerical 
integration  process  which  seems  the  "influence"  for  a  series  of  sources 
as  opposed  to  a  single  one  has  to  be  applied.  Unfortunately  this  signific¬ 
antly  increases  the  number  of  integrations  involved.  The  combination  of 
the  solution  for  the  fluid  domain,  expressed  using  BEM,  the  shell  discretiza 
tion  using  FEM  will  be  the  topic  of  a  future  progress  report. 

The  present  report  deals  with  the  accurate  confutation  of  the  forces 
on  the  semimerged  shells  and  concludes  that  it  is  necessary  to  properly 
define  the  geometry  of  the  obstruction.  Otherwise  an  uneconomic  number 
of  fluid  boundary  elements  is  required  to  obtain  accurate  solutions. 

Higher  order  elements  as  the  one  advocated  in  this  report  are  also  necessary 
to  maintain  compatibility  of  shell  and  fluid  movements. 
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2.  Governing  Equations 

The  following  equations  correspond  to  two  dimensional,  incompressible 
and  irrotational  flow.  The  linear  wave  theory  potential  for  constant 
depth  and  without  the  obstruction  can  be  written  as  (figure  I) 


4>o(x,y,t) 


iga  cosh[K(y+h)l  .  .  _ 

°  o  u  '  J  ux  -twt 

■  .  — M  -  —  III  -i  ■  ■  -  —  II  I  g  p 

a)  cosh  Kh 


(1) 


where  the  subscript  'o'  refers  to  the  incident  field,  x  is  the  coordinate 
in  the  direction  of  the  wave,  y  the  vertical  coordinate  measured  from  the 
mean  water  level,  g  is  gravity,  a^  is  the  wave  amplitude  of  the  incident 
wave,  (o  and  k  the  wave  frequency  and  number  respectively. 

This  wave  potential  can  be  differentiated  in  the  direction  of  the 
normal  to  the  boundaries  of  the  domain  to  give. 


ga^K  r  cosh[x(y+h)J  sinh Jjc(y+h)] 

a)  px  cosh  tch  1  ny  cosh  Kh 


Q0(z»y*t)  =  qo(x,y)e*lwt  (2) 

where  n^  and  n^  are  the  direction  cosines  of  the  normal  to  the  boundary 
with  respect  to  x  and  y. 

The  potential  of  the  scattered  wave  can  now  be  written  as, 

^(x.y.t)  =  4>s  (x,y)e  iu)t  (3) 

The  total  potential  is 

$  *  $  +  $  (4) 

o  s 

and  have  to  satisfy  the  Laplace  equation  in  the  domain,  i.e.  (figure  2) 
Vz4>(x,y,t)  »  0  in  £1  (5) 


with  the  following  b-c 
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i)  Bottom  condition. 
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ii)  Free  surface  condition 
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iii)  Obstruction  condition 
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(7) 


(8) 


iv)  Radiation  -  Sunanerfeld  type  -  condition 


lim 
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(9) 


The  $  solution  for  the  incident  wave  will  satisfy  these  conditions, 
o 

hence  the  problem  can  be  expressed  in  terms  of  the  scattered  field 
carrying  out  the  time  derivatives.  This  gives  the  following  system 


Vz$s(x,y)  =  0  in  ft 


(10) 


with  the  following  boundary  conditions 


i)  Bottom  condition, 
3d> 

a 
3n 
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ii)  Free  surface  condition 
3$ 


3n 


8  (A)  _ 
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(12) 
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iii)  Obstruction  condition 


3$  3<j> 

3n  3n 


-,o  onf, 


(13) 


iv)  Radiation  condition  (approximate) 


3<t> 

tt—  --  i K<p  *  0  on  r, 
3n  s  3 


(14) 


3.  Weighted  Residual  Statement 

The  above  system  of  equations  can  now  be  rewritten  in  a  weighted 
residual  form,  to  minimize  errors  when  using  an  approximate  method  of 
solution  such  as  finite  or  boundary  elements.  Let  us  consider  an 

•  .  3  it* 

arbitrary  function  <p*  and  its  derivative  q*  =  .  Later  on  we 

will  associate  these  functions  with  the  virtual  increment  type  of 
functions  used  in  finite  elements  or  with  the  fundamental  solutions  of 
boundary  elements.  The  resulting  weighted  residual  statement  can  be 
written 


(V2<J>)$*dfl 


+ 


i  «(>)  (J>*dr 


+ 


'•3n 


~  *)$*dr 


(15) 


Notice  that  the  scattered  potential  field  function  $g  is  now  written 
without  the  subscript  's'  for  simplicity. 

Integrating  by  parts  once  the  terms  on  the  left  hand  side  of  (15) 
results  in. 
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q^  <p*dr 


+  i  |  k w*ar  +  j  ad  w* dr  06) 


This  expression  is  the  starting  point  for  finite  element  solutions 
of  the  wave  diffraction  problem. 

If  we  integrate  (16)  once  more  one  obtains. 


4  ?2<f>*dfl  =  - 


‘J 


+  i  I  «p  <fi*dr  + 


<P<P  *dr 


07) 


Notice  that  V  is  the  total  boundary  i.e.  r  =  rj  +  r2  +  ^3  +  r4* 

Expression  (17)  is  the  starting  point  for  boundary  element  solutions. 
We  need  to  define  <f>*  as  a  fundamental  solution  such  that  the  integral 
in  Q  disappears  and  the  problem  becomes  a  boundary  problem  only. 


4 .  Finite  Element  Formulation 

Let  us  first  consider  the  finite  element  discretization  of  equation 
(16).  This  formulation  will  imply  to  define  some  type  of  interpolation 
functions  over  an  element  (figure  3).  The  <}>  function  can  now  be  written 
as, 

♦  -  M  ♦"  (18) 

where  N  are  the  interpolation  functions  chosen  and$n  the  corresponding 
nodal  unknowns.  The  inverted  form  of  (18)  is 

■  6$  *  N  6 q>n 


(19) 


Notice  that  <j>  is  a  complex  function,  i.e. 


<p  m  Ip  +  iip  (20) 

The  resulting  finite  element  matrices  can  now  be  deduced  as 
indicated  in  the  literature  fl"]  and  the  final  system  of  equations  for 
the  continuum  as  shown  in  figure  4  have  the  following  form, 

|tc-(iK+Y)M}$n  =  Q  (21) 

which  can  now  be  solved. 


5 .  Boundary  Element  Formulation 

In  the  boundary  element  case  the  functions  can  be  associated  with 
a  full  space  Green's  function  such  that. 


V2$*  +  aa  =  0 


(22) 


where  A^  is  a  Dirac  delta  function,  whose  integral  is  equal  to  one  at  the 
point  * i ’  and  zero  everywhere  else. 

The  solution  to  (22)  is 


<f>*  = 


(23) 


Substituting  (23)  into  equation  (17)  gives. 


q  $*dr  + 
o 


■  iK$A*dr  +  —  4><#>*dr 

J  8 


♦  |i*  dr  (24) 
r 


Formula  (24)  implies  an  integral  relationship  between  a  point  'i* 
inside  the  Q  domain  and  the  values  on  the  T  boundary  of  the  domain.  If 
the  point  i  is  taken  to  be  on  the  domain  we  have  [2}, 
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♦i 


q  $$*dr  ♦ 
o 


1  * 


i*$$*dr  + 


Q) 

—  4><j>*dr 

g 


♦  Fdr  <25) 

r 


where  c^  *  $  for  smooth  boundaries  and  for  a  sharp  corner  its  value  is 
proportional  to  the  interior  angle  -  can  also  be  deduced  from  consider¬ 
ation  of  constant  field  states,,  see  reference  [2]  jVj . 

Equation  (25)  can  be  rearranged  as  follows. 


c. 
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ri+r2 


3<fr» 

3n 


dr  +  j  ‘ iK ♦*) 


dr 


+ 


<t>*)  <t>  dr 


<t>*q  dr 
o 


0 


(26) 


We  can  now  propose  a  type  of  element  to  discretize  the  r  surface  of 
the  domain  (see  figure  5).  The  type  of  elements  is  of  the  utmost  importance 
in  the  analysis  as  we  will  see  but  for  simplicity  let  us  consider  that 
the  elements  are  constant.  The  surface  r  can  now  be  discretized  into 
N  elements  each  with  a  surface  and  such  that  (26)  becomes. 


c. 
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I 

N  ,+N. 


r . 
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iic<j>*  dr 
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r . 
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1 

<i>*dr 

i 

r. 

j 


0 


(27) 


Note  that  the  total  number  of  elements  is  N,  +  N„  +  N„  +  N.  =  N  and  that 
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the  ’j1  subscript  refers  to  the  element  number. 


We  can  distinguish  two  types  of  integrals. 


h.  . 
ij 


dr 


(28) 


and 


ij 


J  dr 

r. 

j 


These  integrals  represent  influence  functions  between  element  i  at 

which  the  fundamental  solution  is  applied  and  any  other  element  ’j’ 

under  consideration.  Note  that  for  the  case  of  linear  elements  g. .  HO, 

11 

as  n  and  V  are  orthogonal.  The  value  of  the  h. .  coefficient  is  h. .  = 

11  n 

h^  +  c^ .  These  coefficients  can  now  be  arranged  in  matrix  form,  each 
row  of  the  matrix  corresponding  to  a  new  * i *  point.  Equation  (27)  now 
becomes , 


E» 


.  i 


h2  H3-iKe3 


H  -  —  G  1 

A  o  A  J 


!2 

?3 


^ J 


E?,  3{Q0i> 


(29) 


Notice  that  H  and  G  matrices  are  real  but  $  and  Q  are  complex. 


$  =  •(>  +  i  <p 

,n  ',n  .n 


So.  =  h  +  1  ?. 


(30) 


Here  we  can  write, 
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(31) 


Solution  of  liquation  (.31)  will  give  the  values  of  scattered  potential 
required  to  solve  the  problem. 


b .  Wave  Forces 

One  can  now  compute  the  forces  on  the  obstruction  from  Bernoulli  s 
equation.  The  dynamic  pressure  is 


p(x,y,t) 


-  p  Re 


8$ (x,y , t) 
3t 


where 


<f  =  (^  +  i  ip )  e 


-iwt 


Hence 


(32) 


p  =  -  p  Re 


{-iu)( (p  +  iip)e  } 


(33) 


Hence, 


=  p  w(y>sina)t  -  ijjcoswt) 
p  =  pwR  sin(o>t  -  a) 


with  R  =  /p-  +  ip2  ,  tana  =  ^  . 

The  forces  in  the  horizontal  and  vertical  directions  can  now  be 
obtained  integrating  the  pressures,  i.e. 


F  =  - 


p  ( -n) dT 


(35) 


n  is  the  vector  normal  to  the  obstruction;  n  =  (sine,  -cosP)  with 
defined  in  figure  6. 

r  r  r  r 


?  =  J  p(sin6,  -cos0)dT  =  |  j  pdy,  -  j  pdx  j  (36) 
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I 
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lienee  in  discretized  form. 
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dy 
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dx 

(37) 

Nl 
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j  =  i 

p.  Ay 

J  u 

-  ~  I 

j  =  l 

p .  Ax. 

J  J 

We  could  also 


find  moments  about  a  given  part  xoyQ* 


This  gives. 


M  = 
o 


-p  [(x-xo)ny  -  (y-7o)nx]}dr 


(38) 


7 .  Numerical  Results 

The  forces  on  cylindrical  obstructions  of  the  type  shown  in  figure  I 
were  studied  for  two  different  depth  to  radius  ratio,  i.e.  h/a  “  3.0 
and  h/a  -  5.0  and  for  a  range  of  wave  numbers.  Results  obtained  using 
constant  boundary  elements  were  compared  against  linear  finite  elements 
and  the  results  published  by  Chakrabarti  [24j . 

At  first  the  boundary  element  results  were  obtained  by  fitting 
tiic  elements  trying  to  follow  the  curved  surface  of  the  cylinder.  The 
results  obtained  in  this  way  for  a  mesh  of  lb  elements  and  the  cylinder 
were  very  poor  and  are  not  shown  in  the  figures.  Further  investigation 
pointed  out  that  the  reason  for  this  disagreement  was  the  discontinuity 
<>f  the  incident  q  function  between  elements,  specially  for  the  elements 
at  the  bottom.  A  fine  mesh  -  of  200  elements  -  would  be  required  to 

obtain  accurate  results  using  constant  elements.  This  mesh  density  is 
unsuitable  for  three  dimensional  applications  and  points  out  the  need  of 
using  higher  order  elements  which  avoid  discontinuities  at  element  corners. 
The  importance  of  these  discontinuities  is  demonstrated  by  the  reasonable 
accuracy  of  the  BE  solutions  shown  in  figures  7  to  10.  Thee  solutions 
were  obtained  by  using  constant  elements  in  a  step  fashion,  i.e.  the 
elements  joining  at  90°  (or  multiple  of  90°)  angles.  This  grid,  similar 
to  the  one  ou  the  boundary  of  a  finite  difference  grid  gave  accurate 
results  with  a  very  coarse  grid. 


The  finite  element  results  shown  in  figures  7  to  10  were  obtained 
to  compare  against  boundary  elements.  It  is  easy  to  see  that  a  much 
larger  number  of  degrees  of  freedom  would  be  necessary  when  using  finite 
eLi  meats,  specially  for  three  dimensional  problems. 

Figures  11  to  14  compare  linear  boundary  elements  with  Chakrabarti ' s 
results.  These  elements  follow  the  geometry  of  the  half-cylinder  and 
reasonable  agreement  is  obtained  with  46  elements.  It  was  then  decided 
to  study  the  convergence  of  the  results  for  vertical  and  horizontal  forces 
Figures  15  and  16  show  that  the  results  converge  when  the  number  of 
elements  is  increased.  To  inprove  the  rate  of  convergence  one  needs  to 
<.  m|  Loy  quadratic  or  cubic  elements. 

Conclusions 

Three  important  conclusions  were  obtained  from  this  work,  namely 

i)  The  advantage  of  using  boundary  elements  by  comparison  with 
finite  elements  in  order  to  reduce  the  number  of  unknowns.  For 
three  dimensional  problems  this  reduction  is  in  the  order  of  10 
allowing  for  the  solution  of  problems  which  otherwise  would  be 
impossible  to  obtain. 

ii)  The  need  of  using  boundary  elements  which  can  follow  the 
geometry  of  the  structure  under  consideration.  In  particular  it  is 
.idvisabLe  to  use  at  least  quadratic  elements  for  curved  surfaces 

or  even  cubic  elements  which  insure  continuity  of  slopes.  It  is 
important  to  remember  that  as  in  finite  elements,  the  functions 
for  Lhe  potential  should  be  at  least  of  the  same  degree  as  those 
used  for  the  geometry.  Otherwise  appreciable  errors  can  be  intro¬ 
duced  in  the  formulation  | 3j . 

iii)  Future  work  will  involve  the  development  of  a  cubic  boundary 
element  to  make  the  potential  fully  compatible  with  the  type  of 
cubic  functions  used  for  shells. 
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